**********REPLICATION SYNTAX FOR *****
*At the Nexus of Observational and Experimental Research
*Kam, Cindy D. and Marc J. Trussler
*Political Behavior
*Version: Stata v14
*Date: 18 Nov 2016

*tip: Create the directory "C:/PB Replication" so all datasets go there

*****SIMULATIONS
***SD'S (1 1)
***SD FOR U: 2
cd "C:/PB replication"
program modela
version 14
tempname sim
postfile `sim' _b_cons _se_cons _b_T _se_T r2 corrsq case using modela, replace
quietly {
forvalues i = 1/1000 {
drop _all
set obs 500
drawnorm M Z, means(2 0) sds(1 1) corr(1 .5 1) cstorage(lower)
gen T = 0
expand 2
replace T = 1 in 501/1000
gen M_T=M*T
gen Z_T=Z*T
*Case1a
gen true_y = 5 + 2*T + 0*M + 0*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1
gen y = true_y + u
reg y T 
corr y true_y
post `sim' (_b[_cons]) (_se[_cons]) (_b[T]) (_se[T]) (e(r2)) (r(rho)*r(rho)) (1)
*Case1b
drop y true_y u case
gen true_y = 5 + 2*T + .5*M + 1*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1.5
gen y = true_y + u
reg y T 
corr y true_y
post `sim'(_b[_cons]) (_se[_cons]) (_b[T]) (_se[T]) (e(r2)) (r(rho)*r(rho)) (1.5)
*Case2
drop y true_y u case
gen true_y = 5 + 1*T + .5*M + 0*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2
gen y = true_y + u
reg y T 
corr y true_y
post `sim'(_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (e(r2)) (r(rho)*r(rho)) (2) 
*Case2.5
drop y true_y u case
gen true_y = 5 + 1*T + .5*M + 1*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2.5
gen y = true_y + u
reg y T 
corr y true_y
post `sim'(_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (e(r2)) (r(rho)*r(rho)) (2.5)
*Case3
drop y true_y u case
gen true_y = 5 + 1*T + .5*M + 1*Z + 0.5*M_T + 0.25*Z_T
drawnorm u, sds(2)
gen case=3
gen y = true_y + u
reg y T 
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (e(r2)) (r(rho)*r(rho)) (3)
}
}
postclose `sim'
end

set seed 123456789
modela


//Model b
program modelb
version 14
tempname sim
postfile `sim' _b_cons _se_cons _b_T _se_T _b_M _se_M _b_Z _se_Z r2 corrsq case using modelb, replace
quietly {
forvalues i = 1/1000 {
drop _all
set obs 500
drawnorm M Z, means(2 0) sds(1 1)  corr(1 .5 1) cstorage(lower)
gen T = 0
expand 2
replace T = 1 in 501/1000
gen M_T=M*T
gen Z_T=Z*T
*Case1
gen true_y = 5 + 2*T + 0*M + 0*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1
gen y = true_y + u
reg y T M Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (1)
*Case1b
drop y true_y u case
gen true_y = 5 + 2*T + 0.5*M + 1*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1.5
gen y = true_y + u
reg y T M Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (1.5)
*Case2
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 0*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2
gen y = true_y + u
reg y T M Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (2)
*Case2b
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2.5
gen y = true_y + u
reg y T M Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (2.5)
*Case3
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0.25*Z_T
drawnorm u, sds(2)
gen case=3
gen y = true_y + u
reg y T M Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (3)
}
}
postclose `sim'
end

set seed 123456789
modelb


*Model c
program modelc
version 14
tempname sim
postfile `sim' _b_cons _se_cons _b_T _se_T _b_M _se_M _b_M_T _se_M_T r2 corrsq case using modelc, replace
quietly {
forvalues i = 1/1000 {
drop _all
set obs 500
drawnorm M Z, means(2 0) sds(1 1)  corr(1 .5 1) cstorage(lower)
gen T = 0
expand 2
replace T = 1 in 501/1000
gen M_T=M*T
gen Z_T=Z*T
*Case1
gen true_y = 5 + 2*T + 0*M + 0*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1
gen y = true_y + u
reg y T M M_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (e(r2)) (r(rho)*r(rho)) (1)
*Case1b
drop y true_y u case
gen true_y = 5 + 2*T + 0.5*M + 1*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1.5
gen y = true_y + u
reg y T M M_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (e(r2)) (r(rho)*r(rho)) (1.5)
*Case2
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 0*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2
gen y = true_y + u
reg y T M M_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (e(r2)) (r(rho)*r(rho)) (2)
*Case2b
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2.5
gen y = true_y + u
reg y T M M_T
corr y true_y
post `sim'(_b[_cons]) (_se[_cons])   (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (e(r2)) (r(rho)*r(rho)) (2.5)
*Case3
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0.25*Z_T
drawnorm u, sds(2)
gen case=3
gen y = true_y + u
reg y T M M_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (e(r2)) (r(rho)*r(rho)) (3)
}
}
postclose `sim'
end

set seed 123456789
modelc


*Model d
program modeld
version 14
tempname sim
postfile `sim' _b_cons _se_cons _b_T _se_T _b_M _se_M _b_M_T _se_M_T _b_Z _se_Z r2 corrsq case using modeld, replace
quietly {
forvalues i = 1/1000 {
drop _all
set obs 500
drawnorm M Z, means(2 0) sds(1 1)  corr(1 .5 1) cstorage(lower)
gen T = 0
expand 2
replace T = 1 in 501/1000
gen M_T=M*T
gen Z_T=Z*T
*Case1
gen true_y = 5 + 2*T + 0*M + 0*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1
gen y = true_y + u
reg y T M M_T Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (1)
*Case1b
drop y true_y u case
gen true_y = 5 + 2*T + 0.5*M + 1*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1.5
gen y = true_y + u
reg y T M M_T Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (1.5)
*Case2
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 0*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2
gen y = true_y + u
reg y T M M_T Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (2)
*Case2b
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2.5
gen y = true_y + u
reg y T M M_T Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (2.5)
*Case3
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0.25*Z_T
drawnorm u, sds(2)
gen case=3
gen y = true_y + u
reg y T M M_T Z
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (e(r2)) (r(rho)*r(rho)) (3)
}
}
postclose `sim'
end

set seed 123456789
modeld


*Model e
program modele
version 14
tempname sim
postfile `sim' _b_cons _se_cons _b_T _se_T _b_M _se_M _b_M_T _se_M_T _b_Z _se_Z _b_Z_T _se_Z_T r2 corrsq case using modele, replace
quietly {
forvalues i = 1/1000 {
drop _all
set obs 500
drawnorm M Z, means(2 0) sds(1 1)  corr(1 .5 1) cstorage(lower)
gen T = 0
expand 2
replace T = 1 in 501/1000
gen M_T=M*T
gen Z_T=Z*T
*Case1
gen true_y = 5 + 2*T + 0*M + 0*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1
gen y = true_y + u
reg y T M M_T Z Z_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (_b[Z_T]) (_se[Z_T]) (e(r2)) (r(rho)*r(rho)) (1)
*Case1b
drop y true_y u case
gen true_y = 5 + 2*T + 0.5*M + 1*Z + 0*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=1.5
gen y = true_y + u
reg y T M M_T Z Z_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (_b[Z_T]) (_se[Z_T]) (e(r2)) (r(rho)*r(rho)) (1.5)
*Case2
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 0*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2
gen y = true_y + u
reg y T M M_T Z Z_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (_b[Z_T]) (_se[Z_T]) (e(r2)) (r(rho)*r(rho)) (2)
*Case2b
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0*Z_T
drawnorm u, sds(2)
gen case=2.5
gen y = true_y + u
reg y T M M_T Z Z_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (_b[Z_T]) (_se[Z_T]) (e(r2)) (r(rho)*r(rho)) (2.5)
*Case3
drop y true_y u case
gen true_y = 5 + 1*T + 0.5*M + 1*Z + 0.5*M_T + 0.25*Z_T
drawnorm u, sds(2)
gen case=3
gen y = true_y + u
reg y T M M_T Z Z_T
corr y true_y
post `sim' (_b[_cons]) (_se[_cons])  (_b[T]) (_se[T]) (_b[M]) (_se[M]) (_b[M_T]) (_se[M_T]) (_b[Z]) (_se[Z]) (_b[Z_T]) (_se[Z_T]) (e(r2)) (r(rho)*r(rho)) (3)
}
}
postclose `sim'
end

set seed 123456789
modele


///compiling datasets
cd "C:/PB replication"
use modela, clear
gen model = 1
append using modelb
replace model = 2 if model ==.
append using modelc
replace model = 3 if model ==.
append using modeld
replace model = 4 if model ==.
append using modele
replace model = 5 if model ==.
bysort case model: sum corrsq

gen casemodel = (case+model*.01)*100
tab casemodel

//generating Table 2
tabstat _b_T _se_T _b_M _se_M  _b_Z _se_Z _b_M_T _se_M_T _b_Z_T _se_Z_T _b_cons _se_cons, format(%9.3f) by(casemodel) nototal


//generating Table 3
gen true_T=2
replace true_T=1 if case==2|case==2.5|case==3
gen MSE_T = (_b_T-true_T)^2

gen true_M_T = 0
replace true_M_T = 0.5 if case==2|case==2.5|case==3
gen MSE_M_T = (_b_M_T - true_M_T)^2

gen true_Z_T = 0
replace true_Z_T = .25 if case==3
gen MSE_Z_T = (_b_Z_T-true_Z_T)^2

gen true_M = 0
replace true_M = .5 if case~=1
gen MSE_M = (_b_M-true_M)^2

gen true_Z = 0
replace true_Z = 1 if case==1.5|case==2.5|case==3
gen MSE_Z = (_b_Z-true_Z)^2
 
tabstat MSE_T MSE_M MSE_Z MSE_M_T MSE_Z_T , format(%9.2f) by(casemodel) nototal


//Generating Figures
preserve
keep if case==1
histogram _b_T if model==1, name(case1a_A, replace) title(Model a) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==2, name(case1a_B_T, replace) title(Model b) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==3, name(case1a_C_T, replace) title(Model c) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==4, name(case1a_D_T, replace) title(Model d) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==5, name(case1a_E_T, replace) title(Model e) percent xline(2) xlabel(1(1)3) 

histogram _b_M if model==2, name(case1a_B_M, replace) title(Model b) percent xline(0) xlabel(-1(1)1) 
histogram _b_M if model==3, name(case1a_C_M, replace) title(Model c) percent xline(0) xlabel(-1(1)1) 
histogram _b_M if model==4, name(case1a_D_M, replace) title(Model d) percent xline(0) xlabel(-1(1)1) 
histogram _b_M if model==5, name(case1a_E_M, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 

histogram _b_Z if model==2, name(case1a_B_Z, replace) title(Model b) percent xline(0) xlabel(-1(1)1) 
histogram _b_Z if model==4, name(case1a_D_Z, replace) title(Model d) percent xline(0) xlabel(-1(1)1) 
histogram _b_Z if model==5, name(case1a_E_Z, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 

histogram _b_M_T if model==3, name(case1a_C_M_T, replace) title(Model c) percent xline(0) xlabel(-1(1)1) 
histogram _b_M_T if model==4, name(case1a_D_M_T, replace) title(Model d) percent xline(0) xlabel(-1(1)1) 
histogram _b_M_T if model==5, name(case1a_E_M_T, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 

histogram _b_Z_T if model==5, name(case1a_E_Z_T, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 
graph combine case1a_A case1a_B_T case1a_B_M case1a_B_Z case1a_C_T case1a_C_M case1a_C_M_T case1a_D_T case1a_D_M case1a_D_Z case1a_D_M_T case1a_E_T case1a_E_M case1a_E_Z case1a_E_M_T case1a_E_Z_T, col(5) row(5) holes(2 3 4 5 9 10 13 15 20) ycommon name(combining_case1a, replace) ysize(10) xsize(10) title("Case 1")
graph export Case1a.tif, width(1800) replace
//as expected, including irrelevant variables does very little EXCEPT HERE you can see the increase in s.d.
restore

graph drop _all

preserve
keep if case==1.5
histogram _b_T if model==1, name(case1b_A, replace) title(Model a) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==2, name(case1b_B_T, replace) title(Model b) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==3, name(case1b_C_T, replace) title(Model c) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==4, name(case1b_D_T, replace) title(Model d) percent xline(2) xlabel(1(1)3) 
histogram _b_T if model==5, name(case1b_E_T, replace) title(Model e) percent xline(2) xlabel(1(1)3) 

histogram _b_M if model==2, name(case1b_B_M, replace) title(Model b) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M if model==3, name(case1b_C_M, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M if model==4, name(case1b_D_M, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M if model==5, name(case1b_E_M, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5) 

histogram _b_Z if model==2, name(case1b_B_Z, replace) title(Model b) percent xline(1) xlabel(0(1)2) 
histogram _b_Z if model==4, name(case1b_D_Z, replace) title(Model d) percent xline(1) xlabel(0(1)2) 
histogram _b_Z if model==5, name(case1b_E_Z, replace) title(Model e) percent xline(1) xlabel(0(1)2) 

histogram _b_M_T if model==3, name(case1b_C_M_T, replace) title(Model c) percent xline(0) xlabel(-1(1)1) 
histogram _b_M_T if model==4, name(case1b_D_M_T, replace) title(Model d) percent xline(0) xlabel(-1(1)1) 
histogram _b_M_T if model==5, name(case1b_E_M_T, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 

histogram _b_Z_T if model==5, name(case1b_E_Z_T, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 
graph combine case1b_A case1b_B_T case1b_B_M case1b_B_Z case1b_C_T case1b_C_M case1b_C_M_T case1b_D_T case1b_D_M case1b_D_Z case1b_D_M_T case1b_E_T case1b_E_M case1b_E_Z case1b_E_M_T case1b_E_Z_T, col(5) row(5) holes(2 3 4 5 9 10 13 15 20) ycommon name(combining_case1b, replace) ysize(10) xsize(10) title("Case 1'")
graph export case1b.tif, width(1800) replace
restore

preserve
keep if case==2
histogram _b_T if model==1, name(case2a_A, replace) title(Model a) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==2, name(case2a_B_T, replace) title(Model b) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==3, name(case2a_C_T, replace) title(Model c) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==4, name(case2a_D_T, replace) title(Model d) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==5, name(case2a_E_T, replace) title(Model e) percent xline(1) xlabel(0(1)2) 

histogram _b_M if model==2, name(case2a_B_M, replace) title(Model b) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M if model==3, name(case2a_C_M, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M if model==4, name(case2a_D_M, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M if model==5, name(case2a_E_M, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5) 

histogram _b_Z if model==2, name(case2a_B_Z, replace) title(Model b) percent xline(0) xlabel(-1(1)1) 
histogram _b_Z if model==4, name(case2a_D_Z, replace) title(Model d) percent xline(0) xlabel(-1(1)1) 
histogram _b_Z if model==5, name(case2a_E_Z, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 

histogram _b_M_T if model==3, name(case2a_C_M_T, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M_T if model==4, name(case2a_D_M_T, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M_T if model==5, name(case2a_E_M_T, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5) 

histogram _b_Z_T if model==5, name(case2a_E_Z_T, replace) title(Model e) percent xline(0) xlabel(-1(1)1) 
graph combine case2a_A case2a_B_T case2a_B_M case2a_B_Z case2a_C_T case2a_C_M case2a_C_M_T case2a_D_T case2a_D_M case2a_D_Z case2a_D_M_T case2a_E_T case2a_E_M case2a_E_Z case2a_E_M_T case2a_E_Z_T, col(5) row(5) holes(2 3 4 5 9 10 13 15 20) ycommon name(combining_case2a, replace) ysize(10) xsize(10) title("Case 2")
graph export case2a.tif, width(1800) replace
restore


preserve
keep if case==2.5
histogram _b_T if model==1, name(case2b_A, replace) title(Model a) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==2, name(case2b_B_T, replace) title(Model b) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==3, name(case2b_C_T, replace) title(Model c) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==4, name(case2b_D_T, replace) title(Model d) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==5, name(case2b_E_T, replace) title(Model e) percent xline(1) xlabel(0(1)2) 

histogram _b_M if model==2, name(case2b_B_M, replace) title(Model b) percent xline(0.5) xlabel(-.5(1)1.5)
histogram _b_M if model==3, name(case2b_C_M, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5)
histogram _b_M if model==4, name(case2b_D_M, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5)
histogram _b_M if model==5, name(case2b_E_M, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5)

histogram _b_Z if model==2, name(case2b_B_Z, replace) title(Model b) percent xline(1) xlabel(0(1)2) 
histogram _b_Z if model==4, name(case2b_D_Z, replace) title(Model d) percent xline(1) xlabel(0(1)2)  
histogram _b_Z if model==5, name(case2b_E_Z, replace) title(Model e) percent xline(1) xlabel(0(1)2) 

histogram _b_M_T if model==3, name(case2b_C_M_T, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M_T if model==4, name(case2b_D_M_T, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M_T if model==5, name(case2b_E_M_T, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5) 

histogram _b_Z_T if model==5, name(case2b_E_Z_T, replace) title(Model e) percent xline(0) xlabel(-1(1)1)  

graph combine case2b_A case2b_B_T case2b_B_M case2b_B_Z case2b_C_T case2b_C_M case2b_C_M_T case2b_D_T case2b_D_M case2b_D_Z case2b_D_M_T case2b_E_T case2b_E_M case2b_E_Z case2b_E_M_T case2b_E_Z_T, col(5) row(5) holes(2 3 4 5 9 10 13 15 20) ycommon name(combining_case2b, replace) ysize(10) xsize(10) title("Case 2'")
graph export case2b.tif, width(1800) replace
restore


preserve
keep if case==3
histogram _b_T if model==1, name(case3_A, replace) title(Model a) percent xline(1) xlabel(0(1)2) 
histogram _b_T if model==2, name(case3_B_T, replace) title(Model b) percent xline(1)  xlabel(0(1)2) 
histogram _b_T if model==3, name(case3_C_T, replace) title(Model c) percent xline(1)  xlabel(0(1)2) 
histogram _b_T if model==4, name(case3_D_T, replace) title(Model d) percent xline(1)  xlabel(0(1)2) 
histogram _b_T if model==5, name(case3_E_T, replace) title(Model e) percent xline(1)  xlabel(0(1)2) 

histogram _b_M if model==2, name(case3_B_M, replace) title(Model b) percent xline(0.5) xlabel(-.5(1)1.5)
histogram _b_M if model==3, name(case3_C_M, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5)
histogram _b_M if model==4, name(case3_D_M, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5)
histogram _b_M if model==5, name(case3_E_M, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5)

histogram _b_Z if model==2, name(case3_B_Z, replace) title(Model b) percent xline(1) xlabel(0(1)2) 
histogram _b_Z if model==4, name(case3_D_Z, replace) title(Model d) percent xline(1) xlabel(0(1)2) 
histogram _b_Z if model==5, name(case3_E_Z, replace) title(Model e) percent xline(1) xlabel(0(1)2) 

histogram _b_M_T if model==3, name(case3_C_M_T, replace) title(Model c) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M_T if model==4, name(case3_D_M_T, replace) title(Model d) percent xline(0.5) xlabel(-.5(1)1.5) 
histogram _b_M_T if model==5, name(case3_E_M_T, replace) title(Model e) percent xline(0.5) xlabel(-.5(1)1.5) 

histogram _b_Z_T if model==5, name(case3_E_Z_T, replace) title(Model e) percent xline(.25) xlabel(-.75(1)1.25) 

graph combine case3_A case3_B_T case3_B_M case3_B_Z case3_C_T case3_C_M case3_C_M_T case3_D_T case3_D_M case3_D_Z case3_D_M_T case3_E_T case3_E_M case3_E_Z case3_E_M_T case3_E_Z_T, col(5) row(5) holes(2 3 4 5 9 10 13 15 20) ycommon name(combining_case3, replace) ysize(10) xsize(10) title("Case 3")
graph export case3.tif, width(1800) replace
restore

graph drop _all


